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Using the (2+l)-dimensional viscous hydrodynamic code VISH2+1 [1—3], we present systematic 
studies of the dependence of pion and proton transverse momentum spectra and their elliptic flow 
in 200 A GeV Au+Au collisions on the parameters of the hydrodynamic model (thermalization 
time, initial entropy density distribution, decoupling temperature, equation of state, and specific 
shear viscosity r//s). We identify a tension between the slope of the proton spectra, which (within 
hydrodynamic simulations that assume a constant shear viscosity to entropy density ratio) prefer 
larger n/s values, and the slope of the pT-dependence of charged hadron elliptic flow, which prefers 
smaller values of n/s. Changing other model parameters does not appear to permit dissolution of 
this tension. 

PACS numbers: 25.75.-q, 25.75.Dw, 25.75.Ld, 24.10.Nz 



I. INTRODUCTION 

After experiments at the Relativistic Heavy Ion Col- 
lider (RHIC) [4-7] and their theoretical analysis [8-10] 
established that the quark-gluon plasma (QGP) created 
in ultra-relativistic heavy-ion collisions is strongly cou- 
pled and behaves like an almost ideal fluid ("perfect 
liquid" ) with very small viscosity, interest in the theo- 
retical and phenomenological determination of the QGP 
transport parameters, in particular its specific shear vis- 
cosity n/s (i.e. the ratio between its shear viscosity 77 
and entropy density s), soared (see [11, 12] for recent re- 
views). In principle, it should be possible to extract this 
quantity from heavy-ion collision experiments, by com- 
paring the measured hadron spectra and their azimuthal 
anisotropics (in particular their elliptic flow) with theo- 
retical simulations of the collision dynamics which treat 
the QGP shear viscosity as an adjustable parameter 
[13, 14]. In practice, this is a complex and difficult task 
that requires careful and highly constrained simulations 
of all dynamical stages of the collision that sandwich the 
viscous hydrodynamic expansion of the QGP between 
non-equilibrium phases describing (i) the initial geometry 
and early evolution of the fireball before its thermaliza- 
tion and (ii) the final kinetic hadron rescattering stage 
after its hadronization [3, 15]. 

The present work is a contribution to help prepare the 
path for such a phenomenological extraction of (?7/s)qgp- 
It employs viscous hydrodynamics to describe the fire- 
ball evolution, sidestepping the issues related to early 



'Correspond to shen3mps.ohio-state.edu 
tEmail: heinz9mps.ohio-state.edu 
+ Email: huovinenOth.physik.uni-f rankf urt . de 
§ Email: HSongSLBL.gov 



and late non-equilibrium evolution by replacing the out- 
put from the (hypothetical) early non-equilibrium evo- 
lution model by initial conditions for the hydrodynamic 
stage (to be adjusted post facto to final hadron spec- 
tra and multiplicities in central collisions [9]), and the 
late-stage hadronic rescattering and kinetic freeze-out by 
a sudden transition from viscous fluid to free-streaming 
particles, using the Cooper-Frye algorithm [16] along a 
hypcrsurface of constant temperature Td oc - This gener- 
alizes analogous attempts to describe experimental data 
from 200 A GeV Au+Au collisions at RHIC with ideal 
fluid dynamics [9, 17-23] to the case of viscous fluid 
dynamics. Related work has already been reported in 
[13, 24]; what distinguishes the present study from these 
earlier papers is that we use a state-of-the-art equa- 
tion of state that matches the latest Lattice QCD data 
[25, 26] at high temperatures to a realistic, chemically 
non-equilibrated hadron resonance gas at low tempera- 
tures. The construction of this EoS is described in [27], 
except that we here implement chemical freeze-out of the 
stable hadron yield ratios at T c ^ em = 165 MeV, by impos- 
ing appropriate temperature dependent non-equilibrium 
chemical potentials for each hadron species below T c h e m 
[17, 18, 21, 28, 29]. This ensures that the final hadron 
yield ratios from our simulations are consistent with their 
measured values which indicate chemical equilibrium at 
temperature T chcm « 160 - 170 MeV [6, 30, 31]. 

The purpose of this study is not a detailed viscous hy- 
drodynamic fit to the RHIC data; its goal is rather to 
build intuition for systematic trends and parameter de- 
pendences that will be useful in forthcoming more am- 
bitious fit attempts. One feature that disqualifies the 
present model study from being taken too seriously in 
comparison with the experimental data is our assump- 
tion of a constant (i.e. temperature independent) spe- 
cific entropy rj/s. While rj/s is probably small in the 
QGP phase [13, 14, 32], possibly close to the KSS bound 
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(s)kss = Itt ^ * s ex P ecte d to increase dramati- 

cally in the late dilute hadronic phase [35, 36]. This can 
have important consequences for the evolution of flow in 
relativistic heavy-ion collisions [37] which will be studied 
in a separate paper [38]. 

II. HYDRODYNAMIC EQUATIONS, INITIAL 
AND FINAL CONDITIONS 

In this work, we use viscous hydrodynamics to simulate 
the collision system by solving the second-order Israel- 
Stewart equations as described in Ref. [39] . The energy- 
momentum tensor of the fluid is decomposed as 

2 V " / = eu"u v - (p + n)A^' y +7^" (1) 

where e is the local energy density, p is the thermal equi- 
librium pressure (given by the equation of state p(e), see 
below), is the local flow 4-velocity, A**" =g' iv — u^u" 
is the spatial projector in the local fluid rest frame, IT is 
the bulk viscous pressure (which we set to zero in this 
paper, assuming that effects from bulk viscosity can be 
ignored relative to those caused by shear viscosity [40]), 
and it 111 ' is the traceless and symmetric shear pressure 
tensor satisfying u fJl 'K llv = 0. The equations of motion are 
the hydrodynamic equations 

d^ v = 0, (2) 

where denotes the covariant derivative in curvilinear 
(T,x,y,rj) coordinates (see [39, 41] for details), coupled 
to the Israel- Stewart [39, 42-44] evolution equations for 
the viscous pressure components: 

T v 2 t w \r)T J 

The dot on the left hand side stands for the local comov- 
ing time derivative D = u^d^, rj is the shear viscosity, 
a \a> = ^{n u v) j g £ ne velocity shear tensor (see [39, 41] 
for notation), and tv is the microscopic relaxation time 
that controls the evolution of tt^ v (we take r x = 3^ [39]). 

The equations arc solved numerically in the two trans- 
verse spatial directions and time, using the (2+1)- 
dimensional hydrodynamic code VISH2+1 [1-3], assuming 
boost-invariant longitudinal expansion along the beam 
direction. The net baryon density and heat conductivity 
are set to zero. 

To initialize the hydrodynamic evolution we must spec- 
ify the starting time To at which the system is sufficiently 
close to local thermal equilibrium for viscous hydrody- 
namics to be applicable, initial energy density and veloc- 
ity profiles, and the initial viscous pressure tensor n liV . 
We here consider To as a tunable parameter and vary it 
between 0.2 and 0.8fm/c in order to study how it affects 
the final hadron spectra and elliptic flow. 

For the initial energy density profile we study both 
Glauber [2, 45-47] and color glass condensate (CGC- 
fKLN) initializations [48-51], in the optical limit (i.e. 



without accounting for cvent-by-event fluctuations [52- 
56]). Figure 1 shows a comparison of typical initial en- 
ergy density profiles generated from Glauber and CGC 
initializations. In the Glauber model we assume a mix- 
ture of 85% wounded nucleon and 15% binary collision 
contributions to the entropy production [57]. For the 
CGC model we assume that the entropy density is pro- 
portional to the produced gluon density distribution, 
computed with the publicly available fKLN code [58]. 
In central Au+Au collisions, both profiles are normal- 
ized to the same total entropy (adjusted to reproduce 
the total final charged hadron multiplicity dN^/dy in 
these collisions) and converted to energy density using 
the equation of state s95p-PCE (see next section). With 
this normalization, both initializations correctly describe 
the centrality dependence of dN^/dy for ideal fluid dy- 
namics (i.e. for isentropic expansion). 

In the viscous case, viscous heating produces ad- 
ditional entropy, resulting in larger final multiplicities 
which we must correct for by renormalizing the initial 
entropy density profile in such a way that the final multi- 
plicity is held fixed. We perform this renormalization for 
the 5% most central Au+Au collisions (i.e. at b = 2.33 fm) 
and then keep the resulting normalization constant fixed 
for non-central collisions, i.e. we again assume that the 
models produce the correct dependence of initial entropy 
production on collision geometry. It is known, however, 
that the fractional increase of the final entropy over its 
initial value due to viscous heating depends on the size of 
the collision fireball [39] and is therefore expected to be 
larger in peripheral than central Au+Au collisions. For 
the results presented in this paper, we have checked that 
the centrality dependence of viscous entropy production 
is sufficiently weak so that it does not strongly modify 
the centrality dependence of dN^/dy. 

Figure 1 shows that the energy density profile from the 
CGC initialization has a steeper surface gradient than the 
Glauber profile. This leads to larger radial acceleration 
(i.e. radial flow develops more quickly) and is also in part 
responsible for the larger spatial eccentricity of the CGC 
profiles at nonzero impact parameters when compared to 
the Glauber eccentricities [57, 59]. 

The shear viscous pressure tensor tt^ is initial- 
ized with its Navier-Stokes value tt^ v = Irjcr^ where 
Oq V is the velocity shear tensor at time To, cal- 
culated from the initial Bjorken velocity profile, 

= (u T , u x , u y , vP) = (1, 0, 0, 0). 

For the medium's viscous properties the shear viscos- 
ity r\js is the key parameter. According to perturbative 
and lattice QCD, the temperature dependence of rj/s is 
weak over the range explored in heavy-ion collisions at 
RHIC energies. This suggests use of a constant ratio rj/s. 
In this work, the value of rj/s is tuned from 0.08 to 0.24 
in order to study the effects of shear viscosity on the 
hadron spectra and elliptic flow. The influence of a tem- 
perature dependent rj/s will be explored in a forthcoming 
paper [38]. 

Final state hadron spectra are calculated from the hy- 
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FIG. 1: (Color online) A comparison of initial energy density 
profiles at to = 0.4fm/c for "central" (6 = 2.33fm, bottom) 
and "peripheral" (t>=7.5fm, top) Au+Au collisions from the 
Glauber and CGC-fKLN models. Shown are cuts along the 
a;-axis (right panels) and y-axis (left panels). The two profiles 
are normalized to the same total entropy at b — 2.33 fm, using 
the EoS s95p-PCE to convert energy to entropy density. 



drodynamics output via the Cooper-Frye procedure [16] 



E- 



,d 3 N, 



d 3 p (2tt) 3 



P ■ d 3 a(x) fi(x,p), 



(4) 



where £ is the freeze-out surface with normal vector 
d 3 a )1 (x). We take for £ an isothermal surface; cal- 
culations for different freeze-out temperatures are pre- 
sented in Sec. IV B. After computing the spectra of 
all hadronic resonances included in EoS s95p-PCE from 
Eq. (4), we use the resonance decay program [60, 61] from 
the AZHYDRD package [83] to let the unstable resonances 
decay. The pion and proton spectra shown in this work 
include all decay products from strong decays. 

The distribution function on the freeze-out surface can 
be decomposed as / = / eq + Sf into a local equilibrium 
part 



feq(p,x) 



1 



e p-u(x)/T(x) ± i 



(5) 



and a (small) deviation Sf from local equilibrium due 
to shear viscous effects for which we make the quadratic 
ansatz [62, 63] (for other possibilities see [64]) using 



$f{x,p) = f cq {p,x)(lT.f cq {p,x)) 



p p -p v -n llv {x) 



2T 2 (x) (e(x)+p(x)) 



(0) 



(the upper (lower) sign is for fermions (bosons)) for all 
particle species. Sf is proportional to the shear viscous 
pressure tensor ■k i * u {x) on the freeze-out surface and in- 
creases (in our case) quadratically with the particle mo- 
mentum. 



III. EQUATION OF STATE 

To solve the equations (2) and (3) one has to know the 
equation of state p(e) (EoS) of the medium. In this work 
we compare three different equations of state, in order to 
study how the EoS affects the hadron spectra and elliptic 
flow. Two of them, SM-EOS Q [2] and EOS L [39] are 
well known in the literature; the former implements a 
(slightly smoothed) first order phase transition between 
an ideal massless parton gas and a hadron resonance gas 
(HRG), the second is a rough attempt to match lattice 
QCD (LQCD) data [65] above T c to the HRG in a smooth 
crossover transition, as seen in LQCD. In both cases, the 
system is assumed to be in chemical equilibrium all the 
way down to kinetic freeze-out at temperature Jd ec . 

Our third equation of state, s95p-PCE, also interpo- 
lates between the HRG at low temperature and the lat- 
tice EoS at high temperatures, but the matching proce- 
dure is more sophisticated than the one used to construct 
EOS L, and the lattice EoS is based on the recent results 
by the hotQCD collaboration [25, 26]. Furthermore, be- 
low Tchom = 165 MeV temperature, the EoS is that of a 
chemically frozen HRG. The matching procedure using a 
chemically equilibrated HRG is explained in detail in [27] . 
The procedure for the chemically frozen HRG is identical 
since the chemical freeze-out temperature is below the 
temperature where the interpolated EoS deviates from 
the HRG EoS. 

However, the version of s95p-PCE used here deviates 
slightly from the s95p-PCE-vl EoS shown in Appendix 
C of Ref. [27]. First, we have chosen T chcm = 165 MeV for 
the chemical freeze-out temperature, as fitted to experi- 
mental data using thermal models [6, 30, 31] , and we have 
considered as stable particles those with half-life larger 
than 40 fm/c instead of 10 fm/c. Second, our s95p-PCE 
corresponds to a historically slightly earlier stage of the 
parametrization of the EoS than the final version pub- 
lished in Ref. [27]: The fit to the lattice data was done 
without the T = 630 MeV data point. This causes at 
most 0.4% difference between this version and the final 
version of the EoS. We have checked that such a small 
difference docs not cause observable consequences in the 
fluid-dynamical evolution [84]. 

We have built the EoS of the chemically frozen hadron 
gas using the standard procedure in the literature: Be- 
low T c hem the ratios of stable hadron yields are fixed to 
their chemical equilibrium values at Tchcm by finite non- 
equilibrium chemical potentials fJ,i(T) [17, 18, 21, 28, 29]. 
It is worth noting that the ratios of individual parti- 
cle densities are not conserved. What is conserved are 
the ratios of the total densities of stable particles, "hi, 
where total density means the sum of the actual density 
of species i and the additional density of the same species 
that would arise if all unstable resonances in the system 
were allowed to immediately and irreversibly decay. The 
rapid processes that form and decay resonances through 
strong interactions arc still in equilibrium, and thus the 
resonance populations are in equilibrium with the popu- 
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FIG. 2: (Color online) The three equations of state s95p-PCE, 
SM-EOS Q, and EOS L used in this paper. The lower panel 
shows the squared speed of sound cj = |^ as a function of 
energy density e whereas p(e) is shown in the upper panel. 



lations of their daughter particles (see [17, 28] for a de- 
tailed discussion). Thus the chemically frozen system is 
in a state of partial chemical equilibrium (PCE). 

In practice the chemically frozen EoS is evaluated as- 
suming that the evolution is isentropic and the ratios 
fii/s stay constant. Strictly speaking this is not the case 
in viscous hydrodynamics since dissipation causes an in- 
crease in entropy. However, we have checked that in our 
calculations the viscous entropy production from fluid 
cells with temperatures below T c i lcm = 165 McV is small 
(see also the right panel of Fig. 8 in [39] ) and our EoS is 
a good approximation of the physical EoS. 

An analytic parametrization of s95p-PCE is given in 
Appendix A; the EoS can be obtained in a tabulated 
form at [66], where the particles included in the hadron 
resonance gas are also listed. (We included all resonances 
listed in the summary of the 2004 edition of the Review 
of Particle Physics [67] up to 2 GeV mass. Note that our 
s95p-PCE is called s95p-PCE165-vO at [66], to differen- 
tiate it from other versions of the parametrization.) 

The three equations of state are compared in Fig. 2. 
The upper panel shows the pressure and the lower panel 
the squared speed of sound as a function of e. The spike 
in c 2 s (e) at e ~ 0.5GcV/fm 3 results from the sudden 



breaking of chemical equilibrium at T c h C m = 165 McV. It 
has negligible consequences for the expansion dynamics. 
Fig. 2 shows that s95p-PCE is a much softer EoS than 
SM-EOS Q in the QGP phase above T c but much harder 
in the phase transition region around T c . Contrary to 
SM-EOS Q and EOS L, the rapid crossover transition 
between quarks and hadrons that is realized by nature 
doesn't have a well-defined "softest point" [68] which 
would cause the fireball to spend an extended time pe- 
riod in the critical region. Instead, the speed of sound 
never drops much below its value in the HRG, causing 
the fireball to cool rapidly through the phase transition 
[69]. 



IV. SPECTRA AND ELLIPTIC FLOW 

In this section, we discuss the dependence of the trans- 
verse momentum spectra in central 200 A GeV Au+Au 
collisions (0-5% ccntrality, 6 = 2.33fm) and the elliptic 
flow V2(pt) in semiperipheral collisions (20-30% central- 
ity, b = 7.5 fen) for pions, protons and (for v^) all charged 
hadrons, on the EoS and various input parameters dis- 
cussed in Sees. II and III. We have also checked that 
everything we say below about the central collision spec- 
tra also applies, at the same level of precision, to the 
0-averaged spectra in semiperipheral collisions. 

Since the amount of viscous heating depends on the 
input parameters, we retune for each case the normaliza- 
tion of the initial energy density profile in central col- 
lisions so that the same final tt + multiplicity density 
dN^+fdy is obtained. Its value is adjusted by eye such 
that an optimal fit to the measured pion spectrum is ob- 
tained in the low-pr region, px < 1.5GeV/c. As there 
are slight discrepancies between the published data from 
the STAR and PHENIX Collaborations, and these exper- 
iments give their results in different centrality bins, we 
have decided to concentrate on PHENIX results [70, 71] 
when comparing the theoretical curves with experimen- 
tal data. Since we do not attempt to fit these data but 
use the comparison only to illustrate trends, this proce- 
dure is acceptable. A future serious dynamical model fit 
to the data will require proper accounting for systematic 
uncertainties and discrepancies between the different ex- 
periments. 

Since viscous heating effects arc relatively more impor- 
tant in peripheral than in central collisions, our renor- 
malization to constant multiplicities at 6 = 2.33fm leads 
to slightly different pion multiplicities at larger impact 
parameters. For a given EoS, ensuring the same final 
pion multiplicity is equivalent to ensuring the same fi- 
nal total multiplicity. For different equations of state 
(see Sec. IV D) identical pion multiplicities correspond 
to slightly different total multiplicities. 

In the following we show hadron spectra and elliptic 
flow up to transverse momenta of 3GeV/c. We empha- 
size that this is for illustrative purposes only and does not 
imply that we believe hydrodynamics to provide a valid 
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description up to such large pt- When comparing model 
results with experimental data, we judge the quality of 
agreement by focussing on the region pt < 1.5GeV/c 
for pions and pr < 2.5GcV/c for protons (which is 
where we believe hydrodynamics is a reliable approach 
[72]). Specifically for pions, if the calculated spectra 
drop off more steeply than the measured ones above 
Pt = 1.5 GeV/c, we discount this discrepancy, noting that 
this is the region where the experimental spectra begin 
to change from an exponential to a power-law shape due 
to the onset of hard physics. 



A. ry/s-dependence at fixed t = 0.4fm/c and 
T dec = 140 MeV 

Transverse momentum spectra of pions and protons in 
the most central Au+Au collision are shown in Fig. 3 
and in the upper left panel of Fig. 4. The spectra include 
all strong resonance decays. Here we hold initial and fi- 
nal conditions fixed (except for a renormalization of the 
initial peak energy density to ensure the same final multi- 
plicity in all calculations) and vary the specific shear vis- 
cosity rj/s (see figure captions for details). One sees that 
under these conditions larger 77/s-values result in flatter 
spectra; the effect is particularly strong for protons at 
low pt- The main reason is that larger shear viscosity 
leads to larger radial flow, due to a positive contribution 
from it 11 " to the effective transverse pressure gradients at 
early times [41, 73, 74]. 

Figure 3 identifies, however, a second contribution to 
the viscous hardening of the spectra: for ?7/s = 0.16 and 
0.24 and evolution with s95p-PCE, we find that the vis- 
cous correction due to the non-equilibrium deviation Sf 
of the distribution function on the freeze-out surface, 
Eq. (6), is positive for pt > 0.5— 1 GeV/c [85], thus adding 
to the hardening of the spectra from radial flow [86] . This 
is the same sign for Sf as found in [75] (for a different 
EoS) but opposite to what had been found earlier with 
VISH2+1 for smaller values of rj/s using SM-EOS Q (i.e. 
a first order phase transition) [2]. (For rj/s = 0.08 Fig. 3 
shows a negative Sf correction for pions at large pr, of 
same sign but much smaller magnitude than found ear- 
lier [2] with Glauber initial conditions and SM-EOS Q). 
Our finding confirms the fragility of the sign of Sf that 
was already discussed in [2] [87]. 

We note in passing that the positive Sf at large pt 
found here with s95p-PCE is found to be largest in near- 
central collisions (6~0) where it can even lead to a pos- 
itive ^/-correction to the differential elliptic flow V2(pt)- 
At larger b, the S /-contribution to v^{pt) remains nega- 
tive here (see right panels of Fig. 4), as has been consis- 
tently observed in other work [1, 2, 13, 75, 76]. 

The time-integrated effect of the shear viscous pressure 
on the radial flow and the "instantaneous" effect of the 
viscous correction Sf to the distribution function on the 
freeze-out surface together give the total shear viscous 
correction to the hadron spectra. For rj/s = 0.08 we see in 
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FIG. 3: (Color online) Pion spectra for 200 A GeV Au+Au 
collisions at 0-5% centrality from VISH2+1, compared with 
PHENIX data [71]. Results for two different constant values 
of rj/s (0.08 and 0.24) are shown; strong resonance decays are 
included. Solid and dashed lines show the spectra calculated 
from the full distribution function / = / cq + Sf ( "with Sf" ) 
and from the equilibrium part only ("without Sf"). The hy- 
drodynamic evolution starts at ro = 0.4fm/c with an initial 
CGC energy density profile and ends at Td oc = 140 MeV. The 
EoS is s95p-PCE. 



Fig. 3 that with s95p-PCE the Sf correction to the pion 
spectrum is almost negligible, but the upper left panel in 
Fig. 4 shows that the pion and proton spectra are still 
flatter than for the ideal fluid, reflecting the larger radial 
flow caused by the shear viscous increase of the transverse 
pressure gradients [2]. Thus both the effect of viscosity 
on radial flow and Sf contribute to the flattening of the 
hadron spectra. 

Comparing with the experimental data we find that 
both pion and proton spectra favor a relatively large 
shear viscosity, rj/s = 0.16 ~ 0.24. We caution that this 
conclusion is based on calculations done with constant 
(i.e. temperature independent) rj/s and may be subject 
to revision once one properly accounts for the increase of 
rj/s in the dilute late hadronic stage. 

Proceeding to the elliptic flow, we start with a discus- 
sion of the charged hadron V2 in the lower left panel of 
Fig. 4. Here, larger shear viscosity values are seen to 
lead to a stronger suppression of elliptic flow. The right 
panels in Fig. 4 show that this suppression is again the 
consequence of two additive effects: shear viscosity re- 
duces the buildup of anisotropic collective flow, reflected 
in the equilibrium part / cq of the distribution function 
on the freeze-out surface (dotted lines in Fig. 4), but 
the viscous correction Sf causes an additional suppres- 
sion of V2- For T-indcpcndent rj/s, both suppression ef- 
fects increase monotonically with shear viscosity; how- 
ever, the increase of the ^/-correction with rising rj/s is 
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FIG. 4: Upper left panel: Transverse momentum spectra dN/(dy d 2 px) for pions and protons from VISH2+1 for the 5% most 
central Au+Au collisions (6 = 2.33fm), compared with experimental data from the PHENIX Collaboration [71]. Lower left 
panel: Differential elliptic flow V2{pr) for charged hadrons from Au+Au collisions at 20-30% centrality (6 = 7.5fm), compared 
with PHENIX data [70]. Right panels: V2(pr) for pions (top) and protons (bottom). In these panels we compare the elliptic 
flow computed from the full distribution function / = / eq + Sf (dashed, dot-dashed and double-dot-dashed lines) with the 
contribution from the equilibrium part only (dotted lines, "without Sf"). Lines with different symbols show calculations for 
different constant values of the specific shear viscosity n/s, ranging from (ideal hydro, solid lines) to 0.24 as indicated. All 
strong resonance decays are included; charged hadrons comprise ir , K ± , p, p, E ± , E T , H", H + , fT, and Q + . The EoS, initial 
and final conditions are the same as in Fig. 3. 



weaker than that of the viscous suppression of the collec- 
tive flow anisotropy. The stronger suppression of V2 for 
larger rj/s is thus mostly due to the viscous suppression 
of anisotropic flow. 

Since elliptic flow data for identified pions and protons 
in the particular centrality bin shown in Fig. 4 are not 
yet available, we compare in the lower left panel with 
experimental data for unidentified charged hadrons. This 
plot suggests that, even for CGC initial conditions which 
produce more eccentric fireballs than the Glauber model 
[13, 57], the V2 data suggest a smaller value for rj/s, ?]/s = 
0.08 — 0.16, than obtained from the pr-spectra for central 
collisions. 



This tension between the slope of the pr-spectra 
(which tends to favor larger n/s values) and the pt- 
dependence of V2 (which favors smaller values) is generic 
and, as far as we were able to ascertain, cannot be 
resolved with purely hydrodynamic calculations that 
assume constant n/s. A possible solution of this 
problem will likely involve accounting for temperature- 
dependence of n/s and/or the transition to a microscopic 
kinetic description for the late hadronic stage. 
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P T (GeV) P T (GeV) 

FIG. 5: Similar to Fig. 4, but for fixed r;/s = 0.16 and varying decoupling temperature Td oc ranging from 100 to 160 MeV as 
indicated. 



B. Tdec-dependence at fixed To = 0.4fm/c and 
ri/s = 0.16 

In Fig. 5 we explore the sensitivity of spectra and el- 
liptic flow on the value of the decoupling temperature, 
holding all other parameters fixed. For the constant rj/s 
we select rj/s = 0.16 as a compromise between the values 
preferred by the proton spectra and charged hadron V2, 
respectively, in Fig. 4. 

The left upper panel shows that lower freeze-out tem- 
peratures lead to flatter proton spectra. This is a con- 
sequence of additional radial flow built up during the 
extra time the fireball needs to cool down to lower Td ec - 
As is well known [77, 78], the heavier protons receive 
a larger push to higher px from radial flow than the 
lighter pions. Indeed, Fig. 4 shows that the pion spec- 
tra become steeper as Td oc is lowered [79]. Since pi- 
ons are almost massless on the scale of measured trans- 
verse momenta, the inverse slope of their p^-spectrum 



can be approximated by the relativistic blucshift formula 
[77, 78] r s io pc = Tdccy / il^) where (v±) is the average 
radial flow at Td oc - For pions, the steepening effects on 
their spectrum from decreasing Td ec overwhelm the flat- 
tening effects resulting from the associated increase of 
(v±), causing a net softening of the pion spectra for lower 
freeze-out temperatures. 

^From the lower left panel of Fig. 5 one sees that 
lower decoupling temperatures lead to larger elliptic flow 
V2(pt) for charged hadrons. To fully understand this 
systcmatics it is worth comparing charged hadrons to 
the pr-spectra and V2{pt) of pions (upper left and right 
panels, respectively) which dominate the charged hadron 
yield. The observed tendency reflects a combination of 
three effects: 

(i) Since the pr-spectrum of pions (which dominate the 
charged hadrons) gets steeper, even the same hydrody- 
namic momentum anisotropy would lead to a larger slope 
of V2(j>t), to compensate for the lower yield at high px- 
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(ii) Since the fireball hasn't lost all of its eccentricity 
by the time the QGP converts to hadrons [2], additional 
momentum anisotropy is generated during the hadronic 
stage. Lower decoupling temperatures give the system 
time to develop more momentum anisotropy, leading to 
a larger V2 ■ If the pr-spectrum stays unchanged or gets 
steeper (as is the case for pions in Fig. 5), a larger v-i 
must lead to a larger v<i{$>t)- The combination of effects 

(i) and (ii) is seen in the dotted lines in the upper right 
panel, which reflect the hydrodynamic flow anisotropy at 
decoupling, undistorted by viscous corrections Sf to the 
local equilibrium distributions at freeze-out. The effect 

(ii) decreases with increasing rj/s in the hadronic phase 
(not shown here), so the combined effect may be weaker 
than seen in Fig. 5 if viscous hydrodynamics is replaced 
by a microscopic hadron cascade such as UrQMD in the 
hadronic phase. 

(iii) The (negative) viscous corrections from Sf to i>2 
are smaller at lower temperatures, due to the general 
decrease of the viscous pressure components [2]. This 
contributes the largest fraction of the observed increase 
of (pr) with decreasing Tdec, especially at large pr- 

Combining the information from the two left panels in 
Fig. 5 we conclude that both the proton spectra in central 
collisions and charged hadron V2{pt) in peripheral colli- 
sions favor decoupling temperatures near the lower end of 
the window studied here (i.e. Tdoc = 100 MeV works bet- 
ter than Tdec = 140 MeV). The pion spectra are affected 
by variations of Td ec mostly at px ^ 1 — 1.5 GeV/c where 
they fall increasingly below the experimental data as we 
lower Tdec- However, this is also the region where the hy- 
drodynamic description of the pion spectra is known to 
begin to break down [72], due to the gradual transition 
from soft to hard physics which causes the pion spectrum 
to change from an exponential to a power-law shape. Fo- 
cusing therefore on the region px < 1.5(2.5) GeV/c for pi- 
ons (protons), we conclude that a purely hydrodynamic 
description of the experimental data favors freeze-out 
temperatures near 100 MeV. 

The right panels of Fig. 5 show how Td oc affects the 
elliptic flow of different identified hadrons. Charged 
hadrons mostly reflect the behavior of the dominating 
pions whose i>2(pt) increases with decreasing freeze-out 
temperature. But protons behave differently: At low 
Pt < 1 GeV, their elliptic flow decreases with decreasing 
decoupling temperature, while at high px it increases 
with decreasing Td oc - The latter feature reflects the in- 
creasing hydrodynamic momentum anisotropy and de- 
creasing magnitude of the Sf correction, just like it is 
reflected in the pion and charged hadron v%. The de- 
crease of proton vi at low pt, on the other hand, is a 
consequence of having larger radial flow at lower Tdec 
which pushes the protons to larger pr- So, rather than 
thinking of this effect as a decrease of proton vi at fixed 
Pt, we should think of it as shifting the elliptic flow to 
larger pt- 



C. To-dependence at fixed T7/s = 0.16 and 
T dec = 140 MeV 

The upper left panel of Fig. 6 shows that the pion 
and proton spectra react similarly to a change of the 
starting time tq of the hydrodynamic evolution: smaller 
To values lead to more high-px particles, reflecting more 
radial flow. Starting hydrodynamics earlier allows it to 
generate radial flow earlier, and even though this also 
causes the fireball to cool down to Td ec sooner and freeze 
out earlier, the net effect is still a slight increase of the 
average radial flow at freeze-out. 

For soft momenta px < 1.5GcV/c, the effect of tq on 
V2{pt) is negligible. This is true even for protons, show- 
ing that the increase of radial flow with decreasing To is 
a small effect and not enough to visibly push the pro- 
ton v 2 to larger pt- At higher px, the dependence of 
the charged hadron, pion and proton V2 on tq is non- 
monotonic. The right panels of Fig. 6 show that this 
non-monotonic behaviour is the result of two counter- 
acting tendencies which both depend on To monotoni- 
cally: (i) The elliptic flow computed from the local equi- 
librium part / eq of the distribution function at freeze- 
out increases monotonically with increasing To, reflect- 
ing the longer total fireball lifetime (and thus the longer 
time available to build up momentum anisotropy) when 
the hydrodynamic evolution starts later, (ii) The V2- 
supprcssion resulting from the viscous correction Sf at 
freeze-out also increases monotonically with increasing 
To. We don't have a complete understanding of why 
starting (and thus also ending) the hydrodynamics later 
leads to a larger Sf on the decoupling surface; we suspect 
that since the hydrodynamical flow would eventually set- 
tle into a three-dimensional spherically symmetric Hub- 
ble flow with no shear stress, starting earlier leads to a 
stronger transverse flow, and thus to a flow profile which 
is closer to a spherically symmetric flow at the time of 
decoupling. 



D. EoS dependence at fixed To =0.4fm/c, 
tj/s = 0.16, and T dec = 140 MeV 

In Fig. 7 we study the sensitivity of hadron spectra and 
elliptic flow on the equation of state, holding all other 
hydrodynamic parameters fixed (except for the normal- 
ization of the initial energy density profile which is again 
adjusted to ensure constant final multiplicity in central 
Au+Au collisions). We first note that, due to the differ- 
ent chemical composition at hadron freeze-out, the pro- 
ton yields for EOS L and SM-EOS Q are below those 
of s95p-PCE if we hold the pion multiplicity fixed: In 
s95p-PCE we prohibit protons from annihilating on an- 
tibaryons while such annihilation processes are allowed in 
the other two equations of state which assume hadrons in 
chemical equilibrium. To explore flow effects we should 
concentrate on the shape (i.e. inverse slopes) of the pion 
and proton spectra. We see that EOS L produces the flat- 
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P T (GeV) P T (GeV) 

FIG. 6: Similar to Fig. 4, but for fixed r\js = 0.16 and varying starting time to for the hydrodynamic evolution, ranging from 
0.2 to 0.8fm/c as indicated. 



test spectra, followed by SM-EOS Q, whereas the spectra 
from s95p-PCE are steepest. Since all three curves cor- 
respond to the same (constant) freeze-out temperature 
Idee = 140 MeV, these differences can only arise from dif- 
ferent amounts of radial flow or different Sf corrections 
(i.e. different viscous pressure components n^) along the 
freeze-out surface. To separate these two effects we plot- 
ted the spectra calculated without the Sf correction and 
found the same hierarchy. We conclude that, for fixed 
freeze-out temperature, s95p-PCE produces the weakest 
radial flow averaged over the freeze-out surface, EOS L 
generates the strongest flow, with SM-EOS Q falling in 
between. 

The reasons for s95p-PCE generating less radial flow 
than the other two equations of state are complex and 
subtle. The differences in speed of sound during the evo- 
lution largely cancel out (see Ref. [27]). The key dif- 
ference is that, at a fixed freeze-out temperature, the 
chemically frozen HRG embodied in s95p-PCE has a 
considerably larger energy density (ed cc = 0.301 GcV/fm 3 



at Tdcc = 140 MeV) than the chemically equilibrated 
HRG used in EOS L and SM-EOS Q (which has 
edcc = 0.143 GeV/fm 3 at the same temperature) [17], due 
to the larger-than-equilibrium abundances of of baryon- 
antibaryon pairs and mesons that are prohibited from 
annihilating as the system cools below T c h em - So with 
s95p-PCE the fireball reaches the freeze-out point ear- 
lier, and it has a smaller freeze-out radius. It is this 
latter feature which causes the average radial flow along 
the freeze-out surface to be smaller for s95p-PCE than 
for the other two EoS: when plotting the radial velocity 
profiles along the decoupling surface, we found that all 
profiles are approximately linear functions of the radial 
distance r from the center (qualitatively similar to the 
profiles shown in Fig. 4 of Ref. [80]), and that the profile 
for s95p-PCE has the largest slope. However, the average 
radial flow is smallest because for s95p-PCE the average 
over the freeze-out surface is truncated at a smaller max- 
imal r value. 

The charged hadron, pion and proton elliptic flows 
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V2(pt) show quite large sensitivity to the EoS, espe- 
cially at high pt- (We repeat that the hydrodynamic 
spectra should probably not be trusted beyond pr ~ 
2 — 2.5GcV/c, but plotting them out to 3GcV/c makes 
it easier to see what is going on in the calculation.) But 
we see that most of this sensitivity comes in through 
the 8f correction at freeze-out which is particularly large 
for SM-EOS Q. The reason for this is that the first- 
order phase transition leads to large velocity gradients at 
the QGP-to-mixcd-phase and mixed-phase-to-HRG in- 
terfaces [46] which are largely but not completely washed 
out by viscous effects [2] and leave traces on the decou- 
pling surface. Sf effects are weaker with the smoother 
EOS L than with SM-EOS Q even though EOS L gener- 
ates on average more radial flow. 

To discuss the contribution from collective flow 
anisotropics to pion and proton V2{pt) we focus on the 
dotted lines in the right panels of Fig. 7. We see that, 
while s95p-PCE creates less radial flow, it generates a 



larger flow anisotropy (we checked this by direct com- 
putation), resulting in larger V2{pt) for both pions and 
protons than with the other two equations of state. 
For EOS Q it was found in [17, 18, 21] that if the ki- 
netic freeze-out temperature T^cc is adjusted to repro- 
duce the pT-spectra, the correct implementation of chem- 
ical freeze-out at T c h m in the HRG phase increases the 
mass-splitting between v^{pt) of pions and protons at 
low pt- On the other hand, if the freeze-out temperature 
is kept constant, the mass-splitting at low pt decreases 
[17]. Since we have kept the freeze-out temperature fixed 
in our calculations, we see a similar phenomenon here: 
The elliptic flow mass splitting between pions and pro- 
tons is weaker for the chemically frozen s95p-PCE than 
for the chemically equilibrated EOS L and SM-EOS Q. 
This is a consequence of the weaker radial flow generated 
by s95p-PCE. 
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FIG. 8: Similar to Fig. 4, but for fixed ?7/s = 0.16 and different initial energy density profiles (Glauber vs. CGC) as indicated. 
The elliptic flow from an initial CGC density profile is larger than for the Glauber initialization, due its larger initial eccentricity. 



E. Dependence on the shape of the initial energy 
density profile (CGC vs. Glauber) 

We close with a discussion of the influence of the shape 
of the initial energy density profile on the hadron spec- 
tra and elliptic flow, using the Glauber and CGC-fKLN 
models as examples. For an illustration of these profiles 
see Fig. 1. 

The CGC profile is characterized by slightly steeper 
normalized energy density gradients than the Glauber 
profile. According to the Euler equation for ideal fluids 



this leads to larger radial acceleration. Indeed, the up- 
per left panel of Fig. 8 exhibits slightly natter pion and 
proton spectra for CGC-initialized simulations than for 
Glauber initial conditions. 

The elliptic flow coefficients for charged hadrons, pi- 
ons and protons are all significantly larger for the CGC- 



initialized runs than for Glauber initial conditions. This 
is a direct consequence of the well-known larger eccen- 
tricity of the CGC density profiles [13, 50, 57, 59] which 
drives a larger momentum anisotropy. The effect is quali- 
tatively similar for all hadron species. The small amount 
of added radial flow from the CGC initialization that 
we see in the spectra has very little influence on the pt- 
dependence of the proton V2 when compared to the much 
larger effects coming from the larger source eccentricity. 
The suppression of V2 by viscous Sf corrections at freeze- 
out is similar for CGC and Glauber initial conditions, 
being slightly larger in the CGC case. This is presum- 
ably caused by the slightly larger flow velocities (and flow 
velicity gradients) generated by the CGC profile. 

V. CONCLUSIONS 

We have performed a systematic study of the depen- 
dence of the pion and proton transverse momentum spec- 
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tra and their p^-depcndent elliptic flow on the thermal- 
ization time t , initial energy density profile, equation of 
state, freeze-out temperature and specific shear viscosity 
in (2+l)-dimensional viscous hydrodynamic simulations. 
Assuming a temperature independent shear viscosity to 
entropy ratio and CGC initial conditions for the energy 
density profile, we find that the proton p^-spectra mea- 
sured in 200AGeV central Au+Au collisions at RHIC 
favor 77/s-values between 2 and 3 times the KSS bound 
(§)kss = j~ while the p^-slope of the charged hadron el- 
liptic flow prefers smaller values between 1 and 2 times 
the KSS bound. This tension cannot be resolved by dif- 
ferent choices for the other paramaters whose variation 
we studied. Of course, the 77/s values extracted from 
a comparison with simulations using the less eccentric 
Glauber model for the initial energy density profile are 
smaller, but the comparison with the experimental data 
gets worse (the proton spectra come out steeper) and ten- 
sion between the 77/s values preferred by spectra and tj 2 
gets stronger. Lower freeze-out temperatures improve the 
agreement with the data, in particular with the heavy- 
particle (proton) spectra. We saw very little sensitivity 
to the choice of the termalization time To, but for larger 
values of tq we did not allow for the evolution of pre- 
equilibrium radial and elliptic flow, contrary to what is 
expected to happen in reality. The main reason for not 
doing so was that, at this point, we have no theoretical 
control over this pre-equilibrium flow, and we did not 



want to clutter our study by introducing still further pa- 
rameters. If there is a tendency worth mentioning in the 
context of varying to it is that smaller To-values lead to 
somewhat larger radial flow which helps with the descrip- 
tion of heavy hadron spectra. This may, however, also 
be achievable by starting hydrodynamics later, but with 
non-zero initial transverse flow [23, 81, 82]. 

The main objective of this study was to gain an intu- 
itive understanding what reasonable changes in the key 
parameters of a viscous hydrodynamic simulation will do 
to the final hadron spectra and elliptic flow. By also keep- 
ing an eye on the available experimental data we come 
to the conclusion that a purely hydrodynamic descrip- 
tion of the experimental spectra will probably not work, 
at least not with temperature independent 77/s. Real- 
istic variations of 77/s with temperature are the subject 
of a separate study [38] . Based on that study combined 
with the one presented here we believe that giving up 
on a (viscous) hydrodynamic description of the hadron 
resonance gas stage and replacing it by a more reliable 
microscopic approach is unavoidable for a quantitative 
description of the experimental data. 



Appendix A: Analytic parametrization of EoS 
s95p-PCE [66] 



J 



We used the following analytic parametrization for the equation of state s95p-PCE (energy density e and pressure 
p in GeV/fm 3 , entropy density s in fm -3 , temperature T in GeV): 

1. Pressure: 



0.3299 [exp(0.4346e) - 1] 

1.024-1CT 7 • exp(6.041e) + 0.007273 + 0.14578e 
p(e) = { 0.30195 exp(0.31308e) - 0.256232 

0.332e - 0.3223e°- 4585 - 0.00 3 9 06e- exp(-0.05697e) + 0.1167c- 1 - 233 + 0.1436e- exp(-0.9131e) 
0.3327e- 0.3223e - 4585 - 0.003906e- exp(-0.05697e) 



where ei = 0.5028563305441270GeV/fm 3 , e 2 = 1.62 GeV/fm 3 , e 3 = 1.86 GeV/fm 3 , and e 4 ■ 
GeV/fm 3 . 

2. Entropy density: 



e < e\ 
e\ < e < e 2 
e 2 < e < e 3 
e 3 < e < e4 
e > e4 

(Al) 

9.9878355786273545 



12.2304e 1 16849 
11. 9279c 1 - 15635 
0.0580578+ 11.833c 1 - 16187 
t(e)= \ 18.202e- 62.021814- 4.85479 exp(-2.72407-10- 11 e 4 - 54886 ) 
+65.1272e- 128012 exp(-0.00369624e 1 18735 ) 
18.202e - 63.0218 - 4.85479 exp(-2. 72407-10" 
+65.1272e-°- 128012 exp(-0.00369624e 1 - 18735 ) 



ll e 4. 54886' 

4.75253c- 1 - 

") 



e < e\ 

e\ < e < e 2 

e 2 < e < e 3 

e 3 < e < 64 
e > e4 



(A2) 



where e x = 0.1270769021427449GeV/fm 3 , e 2 = 0.4467079524674040GcV/fm 3 , e 3 = 1.9402832534193788GeV/fm 3 , 
and e 4 = 3.7292474570977285 GeV/fm 3 . 

3. Temperature: 

0.203054e - 30679 
(e+p)/s 



T(e) 



e < 0.5143939846236409 GeV/fm 



e > 0.5143939846236409 GcV/fnr 



(A3) 
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